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Abstract 

A rational expansion of the Fermi density operator is proposed. This 
approach allows to calculate efficiently physical properties of fermionic 
systems at finite temperatures without solving an eigenvalue problem. 
Using N evaluations of the Green's function, the Fermi density operator 
can be approximated, subject to a given precision, in the energy interval 
[— /3, oo] with f3 (X N. The presented method may become especially useful 
for electronic structure calculations involving the calculation of charge 
densities, but may also find other applications in e.g. signal processing 
and numerical linear algebra. 
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1 INTRODUCTION 



Quantum systems are most generally described in terms of their density operator 
p. Once p is known, the expectation values of physical quantities are obtained 
as 

{A) = Tr(pA) (1) 

where A is the associated operator of the quantity under consideration. For 
instance, for the calculation of the charge density, A becomes a projector and 
the charge density is simply given by the diagonal elements of p in the site 
representation. In the following we consider fermionic systems in the grand- 
canonical ensemble where 

p(^,T) = /(^^), (2) 

with H, T and p being the Hamiltonian, the temperature and the chemical 
potential respectively, and 

/(^) = (3) 
1 + e-^ 

being the Fermi function. The Fermi function has been studied extensively, 
and effective approximation schemes for the case of scalar arguments, e.g. the 
Sommerfeld expansion, have been developed 

However, for the calculation of e.g. p, one is faced with / applied to operators. 
For large scale applications, it cannot be be switched into the eigen represen- 
tation of H in order to evaluate Eq. ^, since in general full diagonalization of 
H is practically impossible. Because only polynomial and fractional functions 
of operators can be evaluated, corresponding decompositions of f{x) are highly 
desirable. A recent approach is due to Goedecker |^ who proposed to use sys- 
tematically complex line integrals over the Green's function. In the following a 
fractional expansion is presented which does not depend on the calculation of 
line integrals. It will be shown that physical quantities like the charge density, 
which is at the base of many methods in electronic structure calculations, can 
be obtained effectively without solving an eigenvalue problem, necessitating only 
evaluations of the Green's function at selected points. While the method is a 
priory constructed for finite temperatures, it is also well adapted to approximate 
charge densities at zero temperature since the range of the approximation can 
be arbitrarily extended towards lower temperatures. 



2 The fractional expansion 

It is well known, that the Matsubara expansion of the Fermi function, 

oo 

/(x) = l/2-2^^-^-^-— ^ (4) 
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shows very poor convergence properties when truncated to degree m = N. 
Although being the exact fractional series of f{x), the Matsubara expansion is 
therefore not suited for numerical applications. 
Let us consider the function 

e" 

U{x) := f{x-a)f{-x-a) = — — — — , (5) 

2[cosh(a) + cosh(x)J 

which is depicted in Fig. 1 for a — 20. It is readily seen that for sufficiently 
high a > 0, fa{x + a) will approximate f{x) for all x > —a subject to a given 
precision. We now truncate the series in the denominator, 

e" 

fa{x) K gN{x;a) -, (6) 

2pN{x;a) 

where 

^ x^^ 

PN (x; a) = cosh (a) + ^ — — . (7) 

j=o ^ 

It is readily seen than pn{x; a) has no real zeros. For the fractional expansion 
of Eq. ^, we need all zeros z^, i' ~ I, . . . , N oi q{z) :— pn{x; a)\x2=z (see Q). 
For this purpose we define 

yi = l + cosh(a); y,{z) ^ — — , i = 2,...,N-l. (8) 

[2{i- 

Then, it can be seen that 

y(^) = (yiiz), ■■■,yN{z)f 

satisfies a matrix equation Ay^Zi,) — Zj/y(zy), with the N x N matrix A = 



2 + 2cosh(a) if (i, j) = (1,2); 

2l{2l~l) if ^ {1,1 + I), l = 2,...,N 

"'•^ ~ \ 27V (1 - 27V) if z = 7V; 

else. 



(9) 



One easily shows that the z^ are given by the eigenvalues of the matrix A. It 
is well known that the zeros of a given polynomial can be obtained from an 
eigenvalue problem for a related Hessenberg matrix. Goedeckerj^ already has 
proposed to use this fact for the numerical evaluation of all zeros of a polynomial 
as eigenvalues. The usual scheme corresponds to the implicit choice oiyi — 2*~^; 
here the point is to avoid the explicit use of any factorial by using Eq. [s] leading 
to Eq. |. The Zi, can be obtained as eigenvalues with e.g. QR-rotations in a 
numerically stable way; using standard numerical libraries, TV = 40 still yields 
accurate results, and enhanced precision calculations readily allow for larger 7V. 
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However, A'' = 30 will be already sufficient for many applications as will be 
shown. In the following we stick to even N. For convenience, we also chose TV 
and a so that no duplicate zeros are obtained. From the zeros we obtain the 
2A^-zeros of pN{x;a) as ity'i^. The zeros z^-a do behave well, an example 
is plotted in Fig. 2(i). An important although trivial remark must be made 
on the evaluation of polynomials such as q(z) and q'{z), which is to be done 
numerically using a Horner-like scheme in order to avoid any explicit use of the 
factorial for the reason of limited numerical precision, as e.g. 40! already is a 
number with 39 decimal digits. Denoting 



we now may write down the fractional decomposition 

2N 



fa{x) ~ gN{x\a) = ^ _^ . (10) 



As shown in Fig. 2(ii), the coefficients ^i, also behave well. The approximation 
Eq. |l^ converges rapidly. Choosing e.g. iV = 32 and a — 26, the error in 
approximating fa{x) is less than 10^^ for all real x. We have considered the 
symmetric function /„, since we may now exploit the local symmetry of fa{x) 
about the points x = ±a where fa{x) = 0.5. We can approximate successively 
the Fermi function as sum of shifted functions gN{x] a), 

f{x) w (a; + a; a) + gwix + 3a; a) + ... + gN{x + (2M — l)a; a) 

(11) 

in the range [— (2M — l)a,oo]. This is visualized in Fig. 3 for = 32 and 
a — 26, using 2MN = 192 fractional terms. For x oo, the approximation 



Eq. 11 vanishes like x~^^ compared to exponential decay of the Fermi function, 
resulting naturally in a good approximation. For negative x, the validity range 
of the approximation Eq. ^ may be increased by choosing a higher M, i.e., by 
successively adding shifted realizations of gN{x + (2m — l)a; a). 
We note, that the function f(^a)i(.x) (see Eq. H) represents for sufficiently large 
^ a nearly perfect projector on the subspace x € [—a, a]. The presented ratio- 
nal expansion may therefore also find applications in other fields than physics, 
especially when applied to operators. 



3 Application to operators 

The main interest of the approximation Eq. ^ lies in its generalization as op- 
erator equation, replacing x by some Hamiltonian H. Then, the Fermi density 
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operator 



p{kT) 



1 



(12) 



l + exp(ii^) 



2N M 




H - /i + kT[{2 rn - l)a - x^] 



u—l m— 1 



can be approximated efficiently with 2MN evaluations of the Green's function. 
One may furthermore benefit from the fact that the zeros Xi, as well as the 
corresponding ^^-^a come in quartets x^^x*^, —Xi,, — x* if the Zi, are distinct. 

When applied to operators H, the effect of the approximation Eq. |l^is to cut 
off the contributions of states with eigenvalues smaller than e/ = ix—kT{2M—l)a 
(see also Fig. 3). This has no consequences if the spectrum of H is lower bounded 
with no eigenvalues in this domain. There are certainly applications when this 
effect is wanted, e. g. when considering the contributions of sub-bands separately. 

In the example of Fig. 3, M = 3 may be too small for applications involving 
real metals, since the eigenvalue spectrum is covered down to -3.5 eV only at 
room temperature, and a higher M may be needed. However, M = 3 and N = 
32 already is well adapted for e.g. two-dimensional electron gases in mesoscopic 
systems. Assuming an Fermi level of about 15 meV, the Fermi density operator 
can be approximated quite exactly at temperatures down to 1.5 Kelvin. 

The following simple example demonstrates how the total charge density can 
be obtained without solving an eigenvalue problem. Consider the Hamiltonian 



where the on-site energies have been chosen from a uniform random distri- 
bution Cj s]3,5[, and the hopping amplitudes tj^k have been chosen as 1 for 
J, k = {jxiiy)^ , (kx, ky)'^ being nearest neighbors in the two-dimensional plane. 
This Hamiltonian describes free spinless electrons in a discrete two dimensional 
space in presence of a random impurity potential. Hard wall boundary con- 
ditions have been assumed for a system with A^s = 15 x 15 sites, allowing 
conveniently for direct diagonalization. The chemical potential has been fixed 
between the 25th and 26th smallest eigenvalue of iJ, = (A25 -I- X2e)/2, i.e., the 
system is in contact with a heath bath of constant chemical potential. We are 
interested in the total charge density as function of temperature T , especially 
in the hmit of T ^ 0. 

In our case, the considered system is small enough to calculate the charge 
density rij at site j exactly as 





i=l 



were the u(') arc the normalized eigenvectors, thus allowing for direct compari- 



son. 
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We now approximate the Fermi density operator p according to Eq. 
As noted in the introduction, the charge density is given in this case by the 
diagonal elements of p. The results presented in Fig. 4 have been obtained 
using 2MN = 192 evaluations of the Green's function, thus necessitating the 
solution of linear systems of equations only. It is seen that the charge density is 
indeed very well approximated in the domain where Eq. |ll| approximates Eq. |^ 
as discussed previously. The total charge density at zero temperature (which is, 
of course, given by the number of states with energy smaller than the chemical 
potential, i.e., 25 in the present example), is practically identical to the total 
charge density at low temperatures. 

4 CONCLUSIONS 

A fractional approximation of the Fermi density operator has been proposed 
and the necessary concepts have been presented. This method becomes increas- 
ingly appropriate for higher temperatures where the numerical effort decreases. 
However, its range of convergence can be arbitrarily extended towards lower 
temperatures. It is expected to be useful especially for large scale calculations 
at finite temperatures as e.g. investigations of disordered systems, mesoscopic 
systems and electronic structure calculations in general. 

For valuable suggestions I am grateful to K. Maschke. 
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Figure captions 

Fig. 1: The symmetric function fa{x) for a — 20. 

Fig. 2: Positions of the zeros (i) and the fractional coefficients 7^ (ii) in the 
complex plane for the fractional expansion with A'' = 32 and a = 26. 

Fig. 3: Fermi function f(x) and the fractional expansion with Af = 3, A*' = 32, a = 
26. The dotted lines indicate the M — 3 shifted addents. The error in approximating 
f{x) is less than 10~^ for x > —135. 

Fig. 4: Total charge ntot ~ X/jii ''^i function of 9 := kT/{jj. ~ Amin), where Amin is 
the smallest eigenvalue (see text). The error in approximating ntot using Eq. ^ with 
TV = 32, Of = 18 and M = 3 is less than 10"''. 
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